Spatial variability enhances species fitness in stochastic predator— prey interactions 
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We study the influence of spatially varying reaction rates on a spatial stochastic two-species 
Lotka-Volterra lattice model for predator-prey interactions using two-dimensional Monte Carlo 
simulations. The effects of this quenched randomness on population densities, transient oscillations, 
spatial correlations, and invasion fronts are investigated. We find that spatial variability in the 
predation rate results in more localized activity patches, which in turn causes a remarkable increase 
in the asymptotic population densities of both predators and prey, and accelerated front propagation. 
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Understanding biological diversity has been a central 
issue in ecology [l|, 0, H, 0] • In order to understand the 
coexistence of competing species, several simplified 'toy' 
models for the dynamics of few interacting populations 
such as the paradigmatic Lotka-Volterra predator-prey 
model have been investigated. More recently, the crucial 
role of spatial fluctuations and stochasticity in stabilizing 
such systems has been recognized Indeed, stochastic 
predator-prey models @, that consistently account 
for the internal reaction noise yet reduce to the classi- 
cal coupled Lotka-Volterra differential equations in the 
well-mixed mean-field limit have been found to display a 
remarkable wealth of intriguing features [9|: In contrast 
to the regular nonlinear oscillations of the determinis- 
tic Lotka-Volterra model which always entail a return 
to the initial state, these stochastic spatial models yield 
long-lived, but ultimately decaying erratic population os- 
cillations [nj 11, Tl, 13, Hi, 15 1 . In the absence of spatial 
degrees of freedom, these oscillations can be understood 
through a resonant amplification mechanism for stochas- 
tic fluctuations that drastically extends the transient pe- 
riod before the (finite) system finally reaches the absorb- 
ing stationary state (predator extinction) [16j] . In spa- 
tially extended systems, the mean-field Lotka-Volterra 
reaction-diffusion model is known to support traveling 
wave solutions 17, ID, 13] . In corresponding stochastic 



spatial population models, spreading activity fronts in- 
duce persistent correlations between the prey and preda- 
tor species, and further enhance the amplitude and life 
time of local population oscillations 
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In our studies of different stochastic spatial model vari- 
ants for competing predator-prey populations, we have 
found these intriguing spatio-temporal structures and the 
overall features to be remarkably generic and robust with 
respect to even rather drastic changes of the detailed mi- 
croscopic interaction rules [2(J [22| ■ Yet to render these 
models more realistic and relevant for biological systems, 
one must obviously allow for different fitness of the in- 
dividuals as well as spatial variations in the rates that 
describe the population kinetics. In this letter, we ad- 



dress the latter situation by considering the reaction rates 
to be quenched random variables, drawn from truncated 
Gaussian distributions. This model can be interpreted as 
describing a direct environmental influence on the species 
death and reproduction rates such as, e.g., a local vari- 
ability of available resources. 

By means of individual-based stochastic cellular au- 
tomaton Monte Carlo simulations we find that an in- 
creasing spatial variation of the predation interaction 
or species invasion rate (with fixed mean) enhances the 
steady-state population densities (which we take as a 
measure of the species' fitness) of both predators and prey. 
In contrast, mere variations of the predator death and 
prey reproduction rates have very little effect. While a 
simple mean-field averaging over varying predation rates 
does indeed predict a marked stationary density increase, 
it also grossly overestimates cooperative behavior and 
cannot adequately describe our numerical results. In 
fact, we shall argue that the principal fitness enhance- 
ment mechanism rests in the fact that stronger disor- 
der in the predation rate reduces the size of the local- 
ized regions populated by both species, thus amplifying 
the initial local population oscillations and permitting a 
larger number of activity patches in the asymptotic long- 
time limit. Thus, the fitness enhancement of both species 
through spatial variability, notably in the absence of any 
evolutionary adaption processes, is a consequence of the 
emerging dynamical correlations. Remarkably, we find 
that quenched randomness in the predation rates also 
slightly increases the speed of spreading activity fronts. 

We consider a stochastic Lotka-Volterra model on a 
square lattice (typically with 512 x 512 sites) with pe- 
riodic boundary conditions. Individuals of both particle 
species perform random walks through unbiased nearest- 
neighbor hopping (which occurs with probability one, so 
in effect all rates listed below are to be understood as 
relative to the diffusivity D). We allow multiple, essen- 
tially unrestricted lattice site occupation for particles of 
either or both species (the maximum number per site i 
is capped at n.j < 1000). This eliminates the predator 
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extinction transition present in model variants with re- 
stricted site occupation [1, [3 E(j ■ The 'predator' species 
is subject to spontaneous decay A — > with rate p, in 
contrast with the 'prey' particles that may produce off- 
spring B — > 2B with rate a. When individuals of both 
species meet on any lattice site, a prey is 'eaten' and 
the predators simultaneously reproduce, i.e., we imple- 
ment the predation interaction A + B — > 2A with rate 
A. Our dynamical Monte Carlo simulation proceeds with 
random sequential updates; a Monte Carlo step (MCS) 
is completed once on average each particle in the system 
has been moved and had the chance to react [23| . 

Spatial variability is introduced by drawing the reac- 
tion probabilities for each lattice site from normalized 
Gaussian distributions, truncated at the values and 1, 
with fixed mean (in most cases p, = a = A = 0.5) but 
different standard deviations (7 = 0... 0.9. The reaction 
rates therefore constitute quenched random variables. 

The time evolution of the mean predator density 
PA(t) = (nAi{t)), averaged over 50 Monte Carlo simu- 
lation runs with initially randomly placed particles with 
densities pa(Q) — Pb{0) = 1 is shown in Fig. QJa). The 
absolute value of the Fourier transform (taken over the 
full interval of 500 Monte Carlo steps) of this averaged 
signal, \pa(uj)\, is displayed in Fig. [TJb). Here we have 
used uniform rates p = a = 0.5, while the predation 
rate represents a quenched random variable with mean 
A = 0.5 and standard deviation a\ ranging from to 0.9 
[jrjj ]. For these rates, the prey population density (not 
shown) behaves similarly, with an overall phase shift in 
the transient oscillations (2^|, and both densities reach 
practically identical asymptotic density values, see also 
Fig. [21(a). It is evident that increasing spatial variabil- 
ity markedly amplifies the initial population oscillations 
and reduces the relaxation time towards the steady state. 
Remarkably, both predator and prey densities approach 
larger asymptotic values as o~\ is raised. As shown in 
Fig. [H(a), either species gains a remarkable fitness en- 
hancement by ~ 25% in the investigated a\ range. We 
have also studied spatial variations in the predator death 
rate p and the prey birth rate a, with the other rates 
held uniform. In either case we observe merely a minute 
increase in the few percent range of the asymptotic preda- 
tor and prey densities, not nearly as pronounced as the 
effect of spatially varying predation rates. 

The neutrally stable species coexistence fixed point 
of the classical Lotka-Volterra mean-field rate equations 
gives the stationary predator and prey densities as pa = 
a/X and ps = p/X. Presumably therefore, the fitness 
enhancement of both species stems from those regions 
where the predation rate is significantly lower than the 
average. Before we explore local effects in more detail, 
let us first consider a global average over the truncated 
Gaussian predation rate distributions of these mean-field 
stationary densities. The result is depicted in Fig. [U(a) 
along with the simulation data: The 'naive' averaging 




Time 1 \MCS] 



0.00" 

o.oos 

"e 

^ 0.007 
-S 006 



-a 0.003 
< 

S 0.002 

o 

V 

Q. 

i/l 0.001 
0.000 





- - ax = o o 






-- ox = 0.1 






-- a,, =0.2 






a,. = 0.3 






aj = 0.4 






aj = 0.5 






(J;. = 6 






"X = 0.7 






Oi = 0.8 






<T;. - 0.9 













0.10 0.15 
Frequency Qi 



FIG. 1: (a) Time evolution of the predator density pA(t), 
averaged over 50 Monte Carlo simulation runs on a 512 x 
512 square lattice with initial densities pa(0) = ps(0) = 1, 
predator death rate p = 0.5, prey birth rate a — 0.5, and 
mean predation rate A = 0.5, for different variances a\ as in- 
dicated, (b) Signal Fourier transform | pA (uj) \ . (Color online.) 



procedure indeed yields an increase of both stationary 
population densities; however, it predicts a grossly ex- 
aggerated fitness enhancement owing to the fact that 
mean-field approximations tend to overestimate cooper- 
ative effects [25| . We therefore proceed to investigate the 
prominent role of spatial variations and predator-prey 
correlations in the lattice system. 

As one would expect, increasing disorder broadens the 
peak associated with the transient oscillations in the as- 
sociated Fourier signal, reflecting faster relaxation to- 
wards the asymptotic nonequilibrium stationary state. 
Figure [Bb) clearly reveals the roughly threefold increase 
in amplitude of the stochastic nonlinear population os- 
cillations as a\ is raised from to 0.9. By fitting the 
peak envelopes to a Lorentzian shape (which works well 
except in the pure case with o\ = 0), we extracted the 
characteristic relaxation times t^i^a/b = l/r^/s from 
the full widths at half maximum Ta/b a s function of a\, 
see Fig. Hfb). Note the reduction by a factor ~ 2.5 in 
TreiaxA/s as o\ is increased from 0.1 to 0.6. 

The increasing amplitudes of the initial population 
oscillations suggest that the spatial variability in the 
predation rates tends to cluster both species closer to- 
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FIG. 2: Dependence on the variance a\, measured for uniform 
rates jj, = a = 0.5 of (a) the asymptotic mean population 
densities Pa/b{^ ~* oo), compared with the average over the 
mean-field values (dashed-dotted lines, right-hand y-axis), (b) 
the relaxation time T relax7 ft/ B towards the stationary state, (c) 
the predator/prey correlation lengths I cott a/b ~ predators A: 
full lines, prey B: dashed - and the typical species separation 
distance Ity-p- dashed-dotted line; (d) the front speed Wf r0 nt 

of 

the activity rings, obtained for /i = 0.2 and a = 1.0. 



gether, thus enhancing localized population explosions. 
This interpretation is in fact borne out by measuring the 
steady-state equal-time two-point correlation functions 

C a p(x) 



{n al+x npi) - p a p{3 with a, f) = A,B [27 1. 
After again averaging the data over 50 Monte Carlo sim- 
ulation runs, we have extracted the predator and prey 
correlation lengths I coii a/Bi which essentially measure 
the spatial extent of the population patches, as function 
of ax by least-square fits of Caa(x) and Cbb{x) to expo- 
nentials exp(— |x|/Z cor r) at sufficiently large distance |sc|. 
As depicted in Fig. [2jc) , the predator correlation length 
Icon A decreases by ~ 30% from about 3 to 2.1 lattice 
constants as the disorder variance increases, while I CO ttB 
is reduced by ~ 45% from 2.5 to 1.4. 

Similarly, we infer the typical predator-prey separation 
distance lt yp from the cross-correlation function Cab{x), 
which is negative at short distances, but attains a max- 
imum with positive correlation before tending towards 
as \x\ —* oo [i?]] • Here, we define l typ as the location of 
the maximum of Cab{x). The results as function of the 
standard deviation a\, shown in Fig. Etc), closely follow 
the behavior of the correlation lengths, namely rather 
rapidly decreasing from ~ 3 lattice constants to ~ 1.7 as 
a \ is reduced from 0.1 to 0.4. Thus, when the width of 
the distribution of the spatially varying predation rates A 
becomes larger, the ensuing correlated patches of coexist- 
ing predator-prey populations become more localized in 
regions with low values of A. Consequently, a larger num- 
ber of such patches can be accomodated in the system, 
whereby the long-time population densities increase. The 
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FIG. 3: Propagation speed of radially spreading activity 
fronts in the stochastic Lotka-Volterra model with uniform 
rates p = 0.2 and a — as function of the predation rate A. 
The square-root fit is inspired by the mean-field lower bound. 

stabilizing effect of spatial inhomogcncity has recently 
been elaborated in a two-patch predator-prey model of 
diffusively coupled two-dimensional oscillators 28] . 

The classical two-species Lotka-Volterra reaction- 
diffusion equations, i.e., essentially the mean- field rate 
equations supplemented with diffusive spreading, are 



well-known to support traveling wave solutions [J, Il7l . 
18 . Il9j |. whose minimal front speed can be established 
by standard mathematical tools [29|, [3(J 31 1. Beyond 



the mean-field approximation, however, already in single- 
species systems the incorporation of intrinsic reaction 
noise in the computation of wave front pro pagation ve- 
locities is a rather difficult problem 32l 33. 34 
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and there are very few results available for two-species 
models 
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In the initial stage of our simulations, we observe radi- 
ally spreading circular fronts of prey followed by preda- 
tors [26[. We have therefore set out to numerically mea- 
sure the front propagation speed of spreading rings of 
activity, namely prey invading empty regions followed 
by predators feeding on them, in our two-dimensional 
stochastic Lotka-Volterra model. To this end, we set up 
as initial state a circular patch of B particles, one per 
site, of radius 5 lattice constants, and 10 predators A lo- 
cated on the center site of this patch. In this study, we 
have chosen uniform rates \i = 0.2 and a = 1.0, with spa- 
tially varying predation rate with mean A = 0.5. After 
angular averaging to obtain the radial particle concentra- 
tions, we have determined the invading front location by 
searching for the zero of the first derivative of the radial 
prey density. A linear fit of the data with Monte Carlo 
time yields the front speed Ufrontj which is then averaged 
over typically 50 simulation runs. The change of propa- 
gation speed with the disorder variance ax is plotted in 
Fig. E^d). We find a small but noticeable ~ 1% increase 
of the spreading activity front speed as ax is raised from 
to 0.7, which we interpret as essentially a consequence 
of the larger amplitudes of the more localized popula- 
tion fluctuations caused by the spatial variability of the 
predation rate. Our results for spatially homogeneous 
rates are depicted as function of the predation rate in 
Fig. [31 To avoid problems at small A values due to prey 
population explosions, we chose as initial state a sea of 
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unreproductive B particles (5 per site) and 5 predators 
A located on the center of the grid, with p, = 0.2. The 
data can be fitted reasonably well with a square-root ex- 
pression that is motivated by the known lower bound 
v front > \/^Da(X P — At), where p denotes the prey car- 
rying capacity @, [l7|, El]. Here, 4D A = 1, p = 1000, 
but the dimensionless reaction probability A' p ~ A, so 
indeed the fit constants cq and ci should be of order 1, 
but capture fluctuation-induced renormalizations of the 
mean-field parameters, and the additional offset C2 > 
describes the deviation from the lower bound. 

In conclusion, we have employed Monte Carlo sim- 
ulations to investigate a stochastic two-species Lotka- 
Volterra model subject to quenched disorder in the re- 
action rates on a two-dimensional lattice without oc- 
cupation number restrictions. While randomizing the 
prey birth and predator death rates has little effect, spa- 
tial variability in the predation / species invasion rate 
A markedly enhances the asymptotic densities for both 
predator and prey populations. We provide evidence that 
this remarkable fitness increase is caused by disorder- 
induced modifications in the emerging spatio-temporal 
structures: Upon increasing the width of the random 
rate distribution, the typical length scales of both the 
spatial predator-predator and the prey-prey correlations 
is reduced. This results in more localized patches of ac- 
tivity, presumably in the vicinity of regions where the 
local predation rates are smaller than their mean value. 
The system is thus able to accomodate a larger amount of 
populated regions. We also find that spatial variability in 
the predation rate drastically amplifies the initial popula- 
tion oscillations and markedly reduces the time required 
to reach the steady-state configuration. In contrast, the 
front speed of spreading activity rings from a localized 
center is not very strongly affected by the disorder. Yet 
we do observe that the activity fronts accelerate slightly 
upon increasing the variance of the predation rate. 
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